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We present a formulation of the density-functional theory + Hubbard model (DFT+J7) method 
that is self-consistent over the choice of Hubbard projectors used to define the correlated subspaces. 
In order to overcome the arbitrariness in this choice, we propose the use of non-orthogonal gener- 
alized Wannier functions (NGWFs) as projectors for the DFT+f/ correction. We iteratively refine 
these NGWF projectors and, hence, the DFT+U functional, such that the correlated subspaces 
are fully self-consistent with the DFT+!7 ground-state. We discuss the convergence characteristics 
of this algorithm and compare ground-state properties thus computed with those calculated using 
hydrogenic projectors. Our approach is implemented within, but not restricted to, a linear-scaling 
DFT framework, opening the path to DFT+U calculations on systems of unprecedented size. 
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The physics of localized electrons bound to transition 
metal or Lanthanoid ions is important for understand- 
ing and harnessing the behaviour of complex systems 
such as molecular magnets 1 inorganic catalysts 2 and the 
organometallic molecules that facilitate some of the most 
critical chemical reactions in biochemistry 3 . 

Despite its success at predicting ground-state proper- 
ties of materials, Kohn-Sham density-functional theory 
(DFT) 4 fails to describe the physics of such "correlated- 
electron" systems when local or semi-local exchange- 
correlation (XC) functional are used, often predicting 
results that are not only quantitatively but qualitatively 
inconsistent with experiment 5 . The origin of this ap- 
parent failure has been understood since the work of 
Perdew et al. and is related to the unphysical curva- 
ture of the energy functional with respect to electronic 
occupation number 7,8 inherent to such functionals unless 
a self-interaction correction is employed 9 . 

DFT + Hubbard U (DFT+U) 10 is a simple, compu- 
tationally inexpensive method for improving the descrip- 
tion of on-site Coulomb interactions provided by conven- 
tional XC functionals and, hence, for extending the range 
of applicability of DFT to strongly-correlated materials. 

The principle of DFT+U is to divide the system into 
a delocalized, free electron-like part, the "bath", which 
is well-described by conventional XC functionals, and a 
set of "correlated sites" which is not. The XC functional 
for electrons associated with these sites is then explicitly 
augmented with screened Coulomb interactions, the form 
of which are inspired by the Hubbard- model 11 , together 
with a double-counting term to correct for the component 
already included within the XC functional. 

The correlated sites are defined by a set of 3d 
and/or 4/ atomic-like orbitals, or "Hubbard projec- 
tors" , that are chosen a priori. Projector functions that 
are commonly used include hydrogenic wavefunctions 7 , 
maximally-localized Wannier functions 12 , and LMTO- 
type orbitals 5,10 . This arbitrariness constitutes an unsat- 
isfactory, adjustable parameter in the DFT+U method. 



In this article, we present an approach in which the am- 
biguity in the choice of Hubbard projectors is removed, 
and in which they are determined sclf-consistently with 
respect to the DFT+U ground-state. We first outline 
the theoretical framework of our approach, and present 
results of calculations on ligated iron porphyrin. We ex- 
amine the adequacy of hydrogenic orbitals as Hubbard 
projectors and, in particular, the sensitivity of the results 
to the form of these orbitals. We show that optimized 
non-orthogonal generalized Wannier functions (NGWFs) 
provide an unambiguous and natural choice for Hub- 
bard projectors and we introduce a technique for self- 
consistently delineating the subspaces in which correla- 
tion effects play an important role. 

Our implementation is within the framework of linear- 
scaling DFT, however, the same self-consistent projec- 
tor methodology may be applied to any DFT approach 
that solves for localized Wannier-like functions (either 
directly, or indirectly in a post-processing step using 
an interface to a code such as Wannier90 13 ). Fur- 
thermore, our approach may be readily combined with 
recently-proposed methods to calculate U parameters 
from first-principles 7 ' 14 , facilitating entirely parameter- 
free and self-consistent DFT+U calculations. 

The Hubbard energy correction term in DFT+U can 
be interpreted as a functional that penalizes the unphys- 
ical non-integer occupancy of the spatially localized d— 
or /—orbitals, those that are most prone to the spurious 
self-interaction present in standard DFT XC functionals. 

We use a rotationally-invariant correction term, 
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where U^^ represents the screened Coulomb repulsion 
between electrons of spin a, localized on the correlated 
site /. Eq. (1) is, in effect, a penalty functional for devi- 
ation from idempotency of the projection of the single- 
particle density-matrix onto each correlated subspace. 
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The occupancy matrix in the case of a set of M non- 
orthogonal Hubbard projectors |</4,^), m 6 {1,...,M}, 
localized on site /, is given by 

n C0( CT)m < = Y,ft\^\P£ )m '\^), (2) 

where is a Kohn-Sham eigenstate for spin channel a 
with band index i, crystal momentum k and occupancy 
f£\ and Pi /)m ' = \<p$)(<pW™'\ is the Hubbard pro- 
jection operator. The contravariant dual vectors \(p^ m ) 
are related to the covariant projectors through the site- 
centered overlap matrix , = (ifm^ \<p$ ) which is 
a metric on the correlated subspace C^: \(p^ I ' )m ) = 

\ u} ^\n(I)m'm. Q(I)m'm" q(I) _ ijm' 

Our definition of the occupancy matrix differs to that 
of Refs. 15 and has the following desirable properties: 
the expressions are tensorially correct; the energy and 
resulting potential are rotationally invariant; the result- 
ing potential is Hermitian and localized to the correlated 
site; and the trace of the occupancy matrix gives the 
occupancy of the correlated site 16 . The contravariant 
metric 

Q{i)mm' ig ca i cu i a ted only as an inverse of the co- 
variant overlap matrix 0^ m , , therefore, the duals of the 
Hubbard projectors are also localized to the site. As a re- 
sult, and in contrast with previously proposed approaches 
to DFT+J7 models using non-orthogonal projectors, the 
DFT+J7 potential constructed from the tensorially con- 
sistent energy for a given correlated site remains mani- 
festly local to that site. We note that for the special case 
of an orthogonal set of projectors on each site, the pro- 
jection operator is self-adjoint and the above expressions 
reduce to the DFT+J7 correction of Ref. 7. 

Any set of localized functions may, in principle, be 
used as Hubbard projectors with which to define the oc- 
cupancy matrix. Solutions of appropriate orbital sym- 
metry of the hydrogenic Schrddingcr equation, such as 
atomic-like or linear muffin-tin orbitals, are a common 
choice 5,10 ' 14 . These are generally characterized by an ef- 
fective charge Z that determines their spatial diffuseness. 
For a given value of U, results of DFT+C7 calculations 
with different values chosen for Z will not necessarily 
yield the same ground-state properties 14,17 . Notwith- 
standing, hydrogenic orbitals may be inappropriate in 
cases in which the orbitals of the correlated manifold dif- 
fer significantly from atomic wavefunctions. 

In order to obtain accurate occupancies, a set of projec- 
tors is required which adequately accounts for electronic 
hybridization and which, if possible, is defined unambigu- 
ously for the system under study. Wannier functions, in 
particular maximally-localized Wannier functions (ML- 
WFs) 18 , form just such an accurate minimal basis. They 
have been used with good effect to augment DFT with lo- 
calized many-body interactions 19 , and there is numerical 
evidence to suggest that MLWFs constitute the projector 
set which maximizes the U parameter 12 . 

We work with the single-particle density-matrix, 
which is expressed in separable form 20 p(r, r') — 



4'a(^)K al3 (f>p(r') in terms of a localized basis of NG- 
WFs 21 {0 Q (r)}, related to the Kohn-Sham eigenstates 
by a linear transformation tp^\r) = ^ Q Q (r)M^ <7 ^ a . 
The density kernel K al3 = {<f> a \p\<j)P) is the representa- 
tion of the single-particle density operator p in terms of 
the contravariant duals {0 Q (r)} of the NGWFs, which 
satisfy (</> a \<^) = 6%. The NGWFs are in turn expanded 
in terms of a systematic basis of Fourier-Lagrange, or 
psinc 22 , functions. The size of this basis is determined 
by an energy cutoff, akin to a plane- wave kinetic energy 
cutoff, with respect to which calculations arc converged. 
The DFT energy functional is iteratively minimized with 
respect to both the density kernel and the NGWF ex- 
pansion coefficients. The minimization scheme in the 
ONETEP linear-scaling code is detailed in Refs. 23. 

These NGWFs, therefore, are a readily accessible set of 
localized orbitals which are calculated with linear-scaling 
computational cost. Similarly to MLWFs, NGWF cen- 
tres may be used to calculate polarizabilities 16 . Thus, 
in this framework, it is natural to use a localised subset 
of Wannier functions obtained at the end of a ground- 
state calculation, with appropriate orbital character, as 
Hubbard projectors for defining the DFT+U occupancy 
matrix. NGWFs are adapted to their chemical environ- 
ment, reflecting the balance between the competing ten- 
dencies of electron itinerancy and localization in strongly 
correlated systems and, as a result, provide an accurate 
representation of the occupancy of the correlated site. 

We propose a projector self-consistent scheme whereby 
the Hubbard projectors are determined self-consistently 
by iteratively solving for the Kohn-Sham ground-state 
using the Hubbard projectors defined by NGWFs from 
the DFT+C7 ground-state energy calculation of the pre- 
vious iteration. In this way, the Hubbard projectors con- 
verge to those that are optimally adapted for their own 
DFT+J7 ground-state density. This scheme, as we go on 
to show, rapidly and monotonically converges to an un- 
ambiguously defined DFT+l/ ground-state which, for a 
given U parameter, is of lowest energy. In other words, 
the DFT+L/ energy functional is additionally minimized 
with respect to the set of localized NGWF projectors that 
are, at convergence, self-consistent with the DFT+J7 cal- 
culation from which they are determined. 

We applied our method to iron porphyrin (FeP). Met- 
alloporphyrin systems, such as FeP, play an important 
role in biochemistry. The ability of metalloporphyrins to 
bind simple molecules is of interest, particularly in the 
case of FeP which can have a greater affinity for CO and 
NO than O2, resulting in hindrance of respiration. 

We performed fully converged energy minimization 
on FeP, and its complex with carbon monoxide, using 
the ONETEP code 23 . We used spin-polarized DFT+J7 
within the generalized-gradient (GGA) 24 and pseudopo- 
tcntial 25 approximations. An equivalent plane-wave ki- 
netic energy cutoff of 1000 eV was used with a cubic sim- 
ulation cell of side-length 37 A. The NGWFs were spa- 
tially restricted to atom-centered spheres of radius 5.3 A 
and no density kernel truncation was applied. Since the 
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Hydrogenic projector effective charge 

FIG. 1. (Color online) The interaction energy, positive for an 
unbound ligand, of the CO and FeP moieties (top panel) and 
the magnetic dipole moment projected onto the correlated 
manifold of triplet-state FeP (bottom panel). Both are plot- 
ted at various U as a function of the effective charge Z used 
to define the hydrogenic projectors (solid lines) , while dashed 
lines show those quantities calculated with self-consistent 
NGWF Hubbard projectors. Blue lines indicate the binding 
threshold (top) and the ideal projected moment (bottom). 



principal focus of this study was the dependence of com- 
puted DFT+U ground-state properties on variations in 
the Hubbard projectors for a given U value, optimized 
PBE (U = eV) structures were used. 

Shown in Fig. 1, is the interaction energy between FeP 
and CO as an illustration that the binding affinity be- 
tween moieties in DFT+J7 can be strongly influenced by 
the localization of the Hubbard projectors. As can be 
seen, binding affinity is by no means uniquely defined 
when hydrogenic projectors are used, although this may 
be partly compensated by a projector-dependent first- 
principles 7,14 U parameter. At U — 6 eV it varies from 
approximately 0.04 eV to 0.69 eV over the range of Z con- 
sidered; at U = 4 eV the result is even qualitatively am- 
biguous as a function of Z. Using self-consistent NGWF 
projectors (dashed lines) generally results in energetically 
less favourable ligand binding, demonstrating that, for 
a given value of U, NGWF projectors more effectively 
counteract the spurious tendency of GGA functionals to 
over-bind ligands to FeP 26 . Also shown in Fig. 1, the pro- 
jected magnetic dipole moment of FeP in its ground-state 
varies strongly with the value of Z chosen for hydrogenic 
projectors (solid line), with only a narrow range of Z at 
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Wannier function projector self-consistency iteration 

FIG. 2. (Color online) The difference in total energy E and 
the total energy at projector self-consistency -Escf as a func- 
tion of the projector self-consistency iteration. The procedure 
is initialized (iteration 0) with a set of hydrogenic 3d projec- 
tors to construct the correlated subspace, using the Clementi- 
Raimondi 27 effective charge of Z = 11.17 for iron 3d orbitals. 



U = 6 eV giving values that are close to the expected 
2.0 for optimal projectors. Moreover, a pathological 
inconsistency with experiment emerges in that U values 
of sufficient magnitude to achieve the requisite moment 
(for some Z) bring us into the unphysical regime where 
FeP+CO binding is disfavoured. Conversely, the use of 
self-consistent NGWF projectors (dashed) results in a 
projected magnetic moment which lies within the physi- 
cally reasonable range and is rather insensitive to U. 

Fig. 2 demonstrates the stable convergence of the Hub- 
bard projector self-consistency scheme for FeP+CO at 
different values of U. Each data point represents an in- 
dividual variational total-energy minimization, wherein 
the Hubbard projectors are re-constructed from the op- 
timized ground-state NGWFs from the previous itera- 
tion. The energy decreases rapidly as the projectors are 
refined, converging within a small number of iterations. 
This confirms our understanding that the Wannier are 
optimally adapted for the hybridized character of the 
electronic orbitals, while minimizing the energy. In this 
way, more spatially diffuse self-interaction corrections are 
introduced than with purely atomic orbitals, in a compli- 
mentary manner to such methods as DFT+C7 +V 28 which 
allow more general interaction terms between sites. 

Since we re- use the self-consistent density from the pre- 
vious projector iteration to initialise the following itera- 
tion, much fewer NGWF optimization steps are required 
at each successive projector update step. As demon- 
strated in Fig. 3, this results in an overall computational 
effort for achieving projector self-consistency that is only 
a small overhead compared to the conventional approach. 

In order to achieve meaningful insight into the In- 
dependence of bond formation, it is necessary to allow 
for Hubbard projector update consistent with variations 
in the molecular geometry. We stress that ionic force 
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Wannier function projector self-consistency iteration 

FIG. 3. (Color online) The number of NGWF optimization 
steps needed to converge the total energy for each projec- 
tor self-consistency iteration for FeP+CO. Shown inset is the 
convergence of the correlated subspace, as quantified by its 
3d-orbital character. 



expressions are not complicated by the inclusion of self- 
consistent Hubbard projectors, with no additional terms 
appearing over those in conventional DFT+J7. 

In conclusion, we have proposed and demonstrated a 
method within DFT+t/ for obtaining Hubbard projec- 
tors that are uniquely- denned, optimally adapted to their 
chemical environment, and consistent with the DFT+J7 
ground-state density. Our implementation may be in- 
corporated into any method that either solves directly 
for localized Wannier-like states, or which computes 
such states in a post-processing fashion. If combined 
self-consistently with approaches for calculating U from 
first-principles 7 ' 14 , this work opens up the possibility of 
parameter-free DFT+t/ calculations on large systems. 
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